% 求k
clear,clc;
syms t x(t) v(t) k f(t) m g r;
cd1 = [m,g,r;300,9.8,3];%求k时用到的常数
f(t)=v*k*r^2;
eqn = [m*diff(v,t) == m*g-f];
cond = [v(0)==0];
[vt]=dsolve(eqn,cond);
xt = -1*int(vt,t,[0 t])+500;
x_t_01 = subs(xt,cd1(1,:),cd1(2,:));
x_t_01 = matlabFunction(x_t_01);%得到表二下的x(t)
xdata = [0:3:30];
ydata = [500 470 425 372 317 264 215 160 108 55 1];
[resk,res] = lsqcurvefit(x_t_01,rand(1),xdata,ydata);

%测试k值
% for i = 1:length(xdata)
%     testk = x_t_01(resk,xdata(i));
%     vpa(testk)
% end
save fun xt vt m g r t resk k; 
clear;
load fun;
syms n;
cd2 = [m k g;2000/n resk 9.8];%求解r、n时用到的常数
global funcX funcV des;
funcX = subs(xt,cd2(1,:),cd2(2,:));
funcV = subs(vt,cd2(1,:),cd2(2,:));
des = [n r t];%决策变量

lb = [1 2 0];
ub = [100 4 100];
x0 = rand(3,1); 
fmin_x = fmincon(@cost,x0,[],[],[],[],lb,ub,@nlc);
C = cost(fmin_x);
[rg,rh] = nlc(fmin_x);
resdis = [fmin_x';C rg rh]

% nlcon = @(x) [subs(funcV,des,[x(1),x(2),x(3)]);subs(funcX,des,[x(1),x(2),x(3)])];
%[t2f,t2g] = nlc(x);
nlrhs = [0 0];
nle = [-1 0]; % -1 for <=, 0 for ==, +1 >= 
cl = [-inf 0]';
cu = [0 0]';
xtype = 'ICC';
opts = optiset('solver','BONMIN','display','iter')
Opt = opti('fun',@cost,'nl',@opti_nlc,cl,cu,'bounds',lb,ub,'xtype',xtype,'options',opts,'x0',x0)
[opti_x,fval,ef,info]=solve(Opt,x0)


%根据伞个数、伞规格求总费用
function [C]=cost(x)
C1 = 10000;%单个降落伞伞面费用
[rg,rh] = nlc(x);
if x(2)<=2
    C1 = 65;
elseif 2<x(2)&&x(2)<=2.5
    C1 = 170;
elseif 2.5<x(2)&&x(2)<=3
    C1 = 350;
elseif  3<x(2)&&x(2)<=3.5
    C1 = 660;
elseif  3.5<x(2)&&x(2)<=4
    C1 = 1000;
end
%伞个数*(单个伞面价格+绳子数*单位长绳子价格*单条绳子长+单个伞固定费用)
x(1) = fix(x(1))+1;
C = x(1)*(C1+16*4*x(2)*sqrt(2)+200);
% dis = [x';C rg rh]
end

function [g,h]=nlc(x)   
global funcX funcV des;
g(1) = double(subs(funcV,des,[x(1),x(2),x(3)])-20);
g(2) = double(subs(funcX,des,[x(1),x(2),x(3)])); 
h = [];
end

function [nlc]=opti_nlc(x)   
global funcX funcV des;
x(2) = x(2)/2;
g = double(subs(funcV,des,[x(1),x(2),x(3)])-20);
h = double(subs(funcX,des,[x(1),x(2),x(3)])); 
nlc = [g h];
end

%用函数C1=a*r^bl来拟合r与单个伞面价格曲线
function [ydata] = price_func(x,xdata)
    ydata = x(1)*xdata^x*(2)
end








    
    
    
    



